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ABSTRACT 

Research committed by the Langley Research Center 
through 1995 resulting in the HZETRN code provides the 
current basis for shield design methods according to 
NASA STD-3000 (2005). With this new prominence, the 
database, basic numerical procedures, and algorithms 
are being re-examined with new methods of verification 
and validation being implemented to capture a well 
defined algorithm for engineering design processes to be 
used in this early development phase of the Bush 
initiative. This process provides the methodology to 
transform the 1995 HZETRN research code into the 
2005 HZETRN engineering code to be available for 
these early design processes. In this paper, we will 
review the basic derivations including new corrections to 
the codes to insure improved numerical stability and 
provide benchmarks for code verification. 

INTRODUCTION 

Improved spacecraft shield design requires early entry of 
radiation constraints into the design process to maximize 
performance and minimize costs. As a result, we have 
been investigating computational procedures to allow 
shield analysis starting with preliminary design concepts 
through high-fidelity final design models (Wilson et al. 
2003) 1 . Of particular importance is the need to 
implement probabilistic models to account for design 
uncertainties (Wilson et al. 2004a) 2 in the context of 
optimal design processes (Qualls et al. 2003) 3 . These 
requirements need supporting tools with high 
computational efficiency to enable these design 
processes. Only the HZETRN code has so far been 
identified for this purpose within the NASA STD-3000 
(2005) 4 document. As a result, Wilson et al. (2005a) 5 
have prepared a review of past HZETRN code 
development, verification, and validation. Although there 
has been sporadic research to generalize this code over 


the last ten years, evaluation of code status among 
principal users demonstrated drift in the various code 
versions and databases. As a result of this renewed 
interest in HZETRN for future space systems design, 
there is now a systematic effort to advance, verify, and 
validate these codes in preparation to integrate them into 
engineering design processes. The present paper will 
describe the first few months of this renewed systematic 
effort. 

As NASA’s newly defined development spirals are now 
progressing, there is a need to identify a suitable code 
for the early spiral processes. We have chosen the 
1995 HZETRN version for early spiral development, this 
code has had wide distribution and is largely the code 
used in the prior verification and validation processes 
described by Wilson et al. (2005a) 5 . However, an 
assessment of the 1995 HZETRN code status demands 
the re-evaluation of the computational procedures and 
databases to account for software and database drift 
over the last ten years and the inherent limitations of 
these prior computational procedures. In the present 
paper, we examine the basic derivations to identify 
oversights and improve solution stability and accuracy. 
We use both analytical and Monte Carlo benchmarks to 
both identify problems among various versions and will 
distribute these benchmarks with the code package to 
evaluate future compiler and platform dependent 
problems in addition to future drift for whatever reason. 
This revised 1995 HZETRN code, chronicled as the 
2005 HZETRN, is recommended for standard practice in 
an interim time period until a new version incorporating 
more recent research on methods is defined, verified, 
and validated. 

2005 HZETRN 

The atomic and nuclear processes associated with 
space radiation occur over very short time scales 
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(microseconds) compared with the secular variations of 
the space environment allowing use of a time 
independent master equation represented by a steady- 
state Boltzmann description balancing gains and losses 
(Wilson et al. 1 991 ) 6 of the particle fields (Galactic 
Cosmic Rays, GCR, and Solar Particle Events, SPE) 
interacting with the shield material (including the human 
tissues). The basic procedures for reducing this 
equation to a readily soluble numerical process will be 
reviewed and further developed in the context of the 
marching procedures of the 1995 HZETRN algorithm. 

TRANSPORT THEORY 

The specification of the interior environment within a 
spacecraft and evaluation of the effects on the astronaut 
is at the heart of the space radiation protection problem. 
The Langley Research Center has been developing such 
techniques and an in-depth presentation is given Wilson 
et al. (1 991 ) 6 although considerable progress has been 
made since that publication. The relevant transport 
equations are the coupled linear Boltzmann equations 
for a closed convex domain (Fig. 1) derived on the basis 
of conservation principles for the flux density 
(particles/cm 2 -sr-s-A-MeV) </>.(x, Q, E) for particle type j 
as 

Q»V<t> j (x,Q,E)= XJa jk (Q,Q',E,E') </> k (x,Q',E') dO' dE' 
-a.(E)<j>.(x,Q,E) (1) 

where a.(E) and a. k (Q,Q',E,E') are the shield media 
macroscopic cross sections. The o jk (Q,Q',E,E') repre- 
sent all those processes by which type k particles 
moving in direction Q' with energy E' produce a type j 
particle in direction Q with energy E (including decay 
processes). Note that there may be several reactions 
that produce a particular product, and the appropriate 
cross sections for equation (1) are the inclusive ones. 
Exclusive processes are functions of the particle fields 
and may be included once the particle fields are known. 
Note, we will at times loosely refer to <j>(x, Q, E) as either 
flux or fluence and the usage should be clear from the 
context. The time scale of the processes in equation (1) 
are at most on the order of microseconds while time 
scales of boundary conditions are on the order of 
minutes or longer leaving the resulting interior fields in 
equilibrium with the particles at the boundary. 

The total cross section o.(E) with the medium for each 
particle type is 

Oj (E) = <t« (E) + of (E) + o.r (E) (2) 

where the first term refers to collision with atomic 
electrons, the second term is for elastic scattering on the 
nucleus, and the third term describes nuclear reactions 
where we have ignored the minor nuclear inelastic 
processes (excited single particle states) except for low 
energy neutron collisions. The corresponding differential 
cross sections are similarly ordered. Many atomic 


collisions (~10 6 ) occur in a centimeter of ordinary matter, 
whereas ~10 3 nuclear coulomb elastic collisions occur 
per centimeter, while nuclear scattering and reactive 
collisions are separated by a fraction to many 
centimeters depending on energy and particle type. We 
include in the ofE) term the nuclear decay processes. 
Solution methods first use physical perturbations based 
on the ordering of the cross sections with the frequent 
atomic interactions as the first physical perturbation with 
special methods used for neutrons for which atomic 
cross sections are zero. The first physical perturbation 
to be treated is the highly directed atomic collisions with 
mean free paths on the order of micrometers as 
observed in nuclear emulsion. The usual approximation 
is the continuous slowing down approximation leading to 
well specified range-energy relations as shown in Fig. 2 
but neglects the energy straggling that will be studied 
later in the present treatment (see also Wilson et al. 
2002) 7 . The next term is the highly directed multiple 



Fig. 1 Geometric relations of quantities useful in solving 
equation (1). 



Fig. 2 Range of ions in aluminum. 


coulomb scattering and is usually neglected in many 
models but is of great importance in understanding the 
transport of unidirectional ion beams leading to beam 
divergence and so is treated in detail elsewhere 
(Shavers et al. 1993) 8 . The remaining nuclear reactive 
processes have been given main attention in past code 
developments. We will now consider the formal 
development of the relevant equations for further 
consideration. 

Continuous Slowing Down Approximation 

The collisions with atomic electrons preserve the identity 
of the ion and the differential cross sections are given as 

o.«(Q,Q',E,E) = I n <f‘ jn (B) 8(0 -O') 8 Jk 

xS(E + A(’e n - E') (3) 

where n refers to the atomic/molecular excited states 
with excitation energies s n including the continuum. 
Note, the factor A/ 1 results from the units of E of A MeV 
(equivalent unit of MeV/nucleon). Although the atomic/ 
molecular cross-sections &“ (E') are large (= 10" 16 cm 2 ), 
the energy transfers e n are small (= 1-100 eV) compared 
to the particle energy. The atomic/molecular terms of 
equation (1) may be written as 

Xf<f‘ ik (Q,Q',E,E') <1 >/xO',E) dQ'dE'- <f‘.(E) d/x,Q,E 
= I n (f‘. n (E + A/e n ) 8/x.Q, E + A/ £n ) - cf .(E) <l>.(xO, E) 

= I n {<f‘ jn (E) <p/xO, E) + A/e, d E [&" jn (E) </>.(x,a, E) ]} 
-cf' i (E) 0 i (xO,E) + O(eX) 

= d E [S/E) 8/xO, E)1 + 0(e n 2 ) (4) 

where the stopping power S/E) is given as the sum of 
energy transfers and atomic excitation cross sections as 

S/E) = E n £ n (f‘ jn (E) (5) 

The higher order terms of equation (4) are neglected in 
the continuous slowing down approximation (csda) but 
are discussed elsewhere (Tweed et al. 2005) 9 . Evalu- 
ation of the stopping power by equation (5) is 
deceptively simple in that all of the excited states 
including continuum states of the atomic/molecular 
system need to be known. Furthermore, the projectile 
remains a bare ion except at low energies where the 
projectile ion atomic orbital states begin to resonate with 
the electrons of the media leading to electron capture 
and lowering of the ion charge. These issues are further 
discussed in Wilson et al. (1 991 ) 6 and Tai et al. (1997) 10 . 
Equation (1) can be written in the csda as 

Q»V0/xO,E) - A/d E [Sj(E) <j>.(x,Q, E)] 

= If < 7 . JO, O', E, E') 4> k (x,Q',E') dO' dE' 

-a/E)8/xO,E) ( 6 ) 

where the right-hand side of equation (6) excludes the 
atomic/molecular processes now appearing on the left 
as an energy shifting operator in addition to the usual 
drift term. Neutral particles would have null atomic cross 
sections for which the stopping term of equation (6) does 


not appear. Application of csda in both laboratory and 
space shielding has been wide and the resulting errors 
are discussed elsewhere (Tweed et al. 2005) 9 . Equation 
(6) can be rewritten as an integral equation (Wilson 
1977) 11 as 

<t>/x,Q,E) = {S/E.jP/Ej <t>/r(Q,x),Q,E,) 

+Z/ E Ey dE 'A, P/E ')/™f 4n dE dQ<r/Q,Q', E' E”) 
X<P k (x+[R/E) - Rj(E')]Q,Q',E')}/ S J (E)P/E) ( 7 ) 

where r is the point on the boundary connected to x 
along -Q, E r = R/[p - d + RJ, p is the projection of x onto 
Q (Fig. 1), d is the projection of r onto O, R/E) is the 
distance an ion of type j of energy E will travel before 
losing all of its energy to excitation of atomic electrons, 
and P/E) is the probability a type j ion of energy E will 
have a nuclear reaction in coming to rest in the media. 
The usual range-energy relation is given (Fig. 2) by 

R/E) = fAj dE'/S/E') ( 8 ) 

The nuclear attenuation function is given (Fig. 3) by 

P/E) = exp[- fAj <fj (E') dE'/S/E')] (9) 

where the integral domains in equations (8) and (9) 
extend over the full energy range {0, E}. 



Fig. 3 Probability of nuclear reaction as a function of ion 
type and energy. 

Straightahead Approximation 

The approach to a practical solution of equation (7) is to 
develop a progression of solutions from the simple to the 
complex allowing early implementation of high- 
performance computational procedures and establishing 
a converging sequence of approximations with 
established accuracy criteria and means of verification. 
The lowest order approximation using the straightahead 
approximation was guided by the nucleon transport 
studies of Alsmiller et al. (1968) 12 using the Monte Carlo 
methods in which the differential cross sections were 
approximated as 


tfJQ.a'.E.E') « cf jk (E,E') S(Q-O') (1 0) 

resulting in dose and dose equivalent per unit fluence to 
be within the statistical uncertainty of the Monte Carlo 
result obtained using the fully angle dependent cross 
sections. The relation of angular dependent cross 
sections to spacecraft geometry in space application 
was further examined by Wilson and Kandelwal (1974) 13 
using asymptotic expansions about angular divergence 
parameters demonstrating errors in the straightahead 
approximation to be on the order of the square of the 
ratio of distance of divergence to radius of curvature of 
the shield (a small error in most space systems). 

Equations (6) and (7) were examined for HZE ions using 
the following form for the projectile fragmentation cross 
sections as 

< f ik (Q,0',E,E')= cf jk (E') N, exp{ - [Q^E -Q^E'f/2 £jk 2 } 

( 11 ) 

where o r lk (E') is the cross section for producing fragment 
j from ion k, N, is the normalization constant for the 
exponential function, and e jk is the momentum dispersion 
parameter in the reaction. Substituting the interactive 
form of equation (11) into the integral term (target 
fragments treated separately, Wilson 1977) 11 of the 
Boltzmann equation (6) yields 

if & jk (Q,Q',E,E') <t> k (x,Q',E') dQ'dE' 

= I <f Jk (E0 { <t> k (x,Q ,E) 

+ Ed E (p k (x,Q,E) *J[£jk /(2mE)] 

+ Q»d a <f> k (x,Q ,E) Ejk /(2m E)} (12) 

where the second term on the right hand side of 
equation (12) results from corrections in assuming the 
velocity of the ion is preserved in the interaction (Curtis 
and Wilkinson, 1972, Wilson 1977) 14,11 and the third term 
is error resulting from the straightahead assumption 
(Alsmiller et al. 1968, Wilson 1977) 12,11 . The surprising 
result is that the velocity conserving assumption is 
inferior to the straightahead approximation for the nearly 
isotropic space radiation. Under approximations 
examined in equations (4) and (12) there are great 
simplifications in the Boltzmann equation as given below 

Q*V4>/x,Q,E) - Afd E [S.(E) 4>/x,Q, E) ] 

= I cfJE) <t> k (x,Q,E) - cf/E) 4>/x,Q, E) (13) 

which is strictly applicable to the HZE ions (Z > 2) and 
the light ions and neutrons have additional complications 
arising from the broad energy spectra associated with 
their production although the more favorable straight- 
ahead approximation is useful as indicated in equation 
(12) and as found by Alsmiller et al. (1967, 1968) 12,15 in 
their verification process. The corresponding light ion 
(and neutron) Boltzmann equation is 

Q*V4> ] (x,Q,E) - A/ d E [S .(E) 0/x.Q, E) ] 

= If cf jk (E,E') <p k (x,Q,E') dE'- cf/E) d/x,Q,E) (14) 


where we have made use of the straightahead 
approximation as given by equation (10). Equations (13) 
and (14) have sufficient simplicity to allow an approach 
for both space and laboratory applications. The main 
force of the laboratory applications allow detailed model 
testing of the many atomic/molecular and nuclear 
processes. 

Marching Procedures and HZETRN 

The first version of the HEZTRN code is based on the 
solution of equation (14) as guided by the Monte Carlo 
studies of the straightahead approximation by Alsmiller 
et al. (1968) 12 . We again specialize to solution along a 
ray Q directed along the x-axis for which the Boltzmann 
equation becomes 

M ' Afd E S.(E) + aj 4>/x, E) = ij o jk (E,E') 4>fx,E') dE' 

(15) 

where o jk (E,E') are approximated for nucleons by the 
multiplicities of Ranft (1980) 16 , Bertini et al. (1972) 17 , and 
quasi-elastic contribution as described by Wilson et al. 
(1 991 ) 6 . An immediate problem is the near singular 
nature of the differential operator and transformation 
from energy to residual range coordinates as we did in 
developing the Green’s function greatly relieves this 
problem (Wilson et al. 1 991 ) 6 . Unlike the Green’s 
function development, numerical procedures are 
simplified by introducing only a single residual range 
coordinate for all ions and the residual proton range r is 
used as the common coordinate as 

r = f 0 E dEVS(E') (16) 

and residual range of other particle types is related as 
Vj r k <= r which fails at low energies corresponding to low 
residual range due to electron capture into atomic 
orbitals characteristic to each ion type. The corre- 
sponding transport equation is 

R - Vj d r + ofr)] y/j(x,r) 

= Zj r “ (Vj/v k )s jk ( r ,r') V k (x,r') dr' (17) 

where scaled flux is now (Vjfor neutral particles such as 
neutrons are taken as unity in scaling relations, Wilson 
et al. 1 991 ) 6 

¥j(x,r) - V| S(E) (.(x, E) (18) 

and the scaled differential cross sections are 

s jk (r,r') = S(E)<y jk (E,E') (19) 

Errors in scaling of proton stopping and range 
parameters in arriving at the approximate transport 
equation (17) are compensated in part by solutions of 
equation (17) approaching a low energy equilibrium 
spectrum for ions given by 

v, S(E) 4>.(x, E) => constant (20) 
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where the constant is fixed by the higher ion energy. In 
distinction, the solution to equation (15) for ions has the 
low energy equilibrium spectrum 

A/ 1 Sj(E) 4>.(x, E) => constant (21 ) 

also fixed by the higher energy flux for which the range 
scaling relation Vj rj ~ r has better validity and the two 
constants are nearly equal so that equation (21) has 
improved accuracy over equation (20) at lower energies. 
This fact will require us to alter the flux unsealing 
relations as demanded by equation (21) to maintain 
accuracy at the lower energies. From equations (20) 
and (21), we can understand the simplicity of numerically 
solving equation (17) over numerical solution based on 
equation (15). The solution to equation (17) approaches 
a constant at small residual ranges allowing large 
separations in r grid values with smooth extrapolation to 
zero range whereas solutions of equation (15) vary as 
the nearly singular 1/Sj(E) for which small E grid spacing 
is required leading to slow computational procedures. 
We will test the assumptions in equation (17) and 
unsealing according to relation (21) later in this report. 

The confusion caused by different scaling methods and 
associated coordinates for numerical procedures is 
justified by the simplification of the numerical 
representation of fluence of all particle types over a 
common residual range grid and simplification of the 
numerical procedures leading to high performance 
codes. Still a straightforward finite differencing of 

equation (17) can introduce unstable roots as had 
plagued the thermal transport problem for many years 
(Wilson et al. 1 991 ) 6 . We will use the unconditionally 
stable methods of Wilson et al. (1 991 ) 6 arrived at by 
inverting the differential operator (Wilson et al. 1977, 
1 989,1 991 ) 11,18,6 of equation (7) as 

y/j(x,r) = exp[-Cj(r,x)] i///0,r + Vj x) 

+ zJo x /r + v,x'°°exp[-£i(r,x')](v/vJ s jk ( r +Vj x\ r ') 
x y/ k (x - x',r') dr' dx' (22) 

where the exponential is the integrating factor related to 
attenuation of the j type ions with 

t;j(r,x)=J*o j (r + v j x')dx' (23) 

and is related to equation (9). Equation (22) is a Volterra 
equation and can be solved either as a Neumann series 
or with marching procedures. Note that the inverse 
mapping is taken as 

<j>.(x,E)=A. ¥] (x,r)/S.(E) (24) 

to guarantee the equilibrium solution given as equation 
(21) at low energies away from the boundaries (note, the 
proton stopping power is used in case of unsealing the 
neutron flux). We assume in this treatment that the 
equilibrium constant resulting from equation (22) and 
given in equation (20) differs little from condition (21) for 
which the inverse mapping of equation (24) is most 


accurate. We will be testing these approximations later 
in this report. 

Two tracks are taken in implementing a 
marching procedure for equation (22) depending on 
particle type as demanded by the character of the 
nuclear processes. The problem naturally divides into 
“light ions” which will refer to all ions with atomic mass of 
four or less including neutrons and high charge-energy 
(HZE) ions have atomic mass greater than 4. The 
distinction arises from the energy and angle distributions 
of the double differential cross sections for which the 
HZE ions leaving a projectile fragmentation event have 
velocity nearly equal to that of the projectile as 
approximated by equation (11). Although the light ions 
are assumed to travel in the same direction as the 
projectile (see equation 10), they cover a broad energy 
distribution that cannot be ignored. The marching 
procedure is obtained by first considering equation (22) 
evaluated at x + h where h is the step size. Clearly, 

V fj(x+h,r) = exp[-C(r,h)] y/j(x,r+Vjh) 

+ Zjf L vj/ exp[-C/r,x')](v j/v k )s jk ( r +Vj x\ r ') 
x \j/ k (x+h -x',r') dr' dx' ( 25 ) 

which may be used to develop a marching step from xto 
x + h once a means to approximate the field function 
i j/j(x,r) across the subinterval {x, x+h}. If h is sufficiently 
small such that 

o/r) h« 1 (26) 

then, following lowest order perturbation theory (Wilson 
and Lamkin 1975) 19 , one has 

V r k (x +h - x',r') = exp[-C k (r',h - x')] y/ k [x, r' + v k (h-x')] 

+0(h) (27) 

which may be used to approximate the integral in 

equation (25) giving results for the fields O(tf) as 
required to control the propagated error (Wilson et al. 
1 991 ) 6 . Substituting equation (27) into (25) and 
evaluating the attenuation factors at the interval midpoint 
(mean value theorem) results in 

yz/x+hs) = exp[-C/r,h) ] y/j(x,r+Vj h) 

+ Z k expKj(r,h/2)-Ur',h/2)] 
x J r +y j h/ 2 ° (Vj/vJ F A j k (h,r,r) y/ k (x,r'+v k h/2)]dr' 
+ O(hF) (28) 

where the integrand has been simplified using 

F A jk (h,r,r) = J 0 h s jk ( r +vj x\ r ') dx' 

= F Jk (r + Vj h, r) - F Jk (r,r) (29) 

and 

F jk (r,r) = /„*> O jk (E",E') dE" (30) 

with e(r) the energy associated with proton residual 
range r, and E' = e(r'). Note that if j corresponds to a 
neutral particle such as the neutron (J = n) then the 
above expressions are evaluated in the limit as v s 
approaches zero in the range scaling relations resulting 
in the following (whereas the flux scaling factor for 
neutrons assumes v„ = 1) 



¥n(x + h,r) = exp[-o n (r) h] y/ n (x,r) 

+ Z k * n exp[-o n (r) h/2 -£ k (r',h/2)] h 
xJr (l/v k )s nk (r,r') Vk(x, r'+v k h/2) dr' 

+ exp[-a n (r) h/2 -a n (r')h/2] h f r °° sj r , r ') y/ n (x, r') dr' ( 31 ) 
and similarly for the neutral k term (k - n) when the j 
particle is charged 

¥(x+h,r) = exp[-^(r,h)]y/j(x,r+Vih) 

+ Z k , n exp[-Cj(r,h/2) -Ur',h&)]/ r "WJF*jfary+Vjh&) 
xy/Jx, r'+(vj+v k )h/2)] dr' + exp[-^(r,h/2) -a n (r') h/2] 
xj r “ Vj F A jn (h,r, r '+ v f h/2) y/ n (x, r'+v ) h/2) dr' (32) 

where v n in the flux scaling relation (24) is taken as unity. 
Equations (31) and (32) are solved on an equally space 
x-grid Ax - h apart and a logarithmic spaced r-grid on 
two subintervals. The remaining integrals in these 
equations are approximated by (Wilson et al 1 991 ) 6 

lk~K(r n ,r ) Wj(x,r') dr' »Z., =k K[r„,(ri+r l+1 )/2] J/ +I y/jfx.r') dd (33) 

where °° denotes a chosen upper limit tailored to the 
specific boundary condition. Note that the matrix of K- 
values can be evaluated once on the r-grid and stored 
for subsequent steps providing high computational 
efficiency. Equations (31) and (32) provide the basis of 
the light ion transport of both the HZETRN and the 
BRYNTRN codes. The HZE ion projectile (Aj > 4) 
coupling to the light fragments is contained in equations 
(28) to (32). 

The HZE fragments are produced with nearly the same 
velocity as the projectile ion as expressed in equation 
(13) and results in the simplified Boltzmann equation as 

M - Af^S/E) + o .(E)] 4>/x, E) = I k o jk (E) p k (x,E) (34) 

for which the scaled equations result in contributions 
from all HZE ions (with A k > 4 ) as 

Vi(x,r) = exp[-C/r,x) ] y//0,r + Vj x) 

+ zJq exp[-£j(r,x')](v/vj 
x cyrfVj x') y/ k (x - x',r + Vj x') dx' (35) 

The corresponding marching equation (Wilson et al. 
1 991 , Shinn et al. 1 992) 6,20 is given as 

¥j( x +h,r) - exp[-Cj(r,h)] y/j(x,r+Vjh) 

+ Z k f 0 h exp[-C/r,x')](v/v k ) 
x a jk (r+VjX') y/ k (x+h-x', r+ vjx') dx' (36) 
for which the integrand can be approximated for 
sufficiently small h using 

¥k(x +h - x',r+VjX') - exp[-C k (r +Vj x',h - x')] 

x y/Jx, r+VjX' + v k (h-x')] +0(h) (37) 

allowing the following simplification 

y/j(x+h,r) - exp[-Cj(r,h)] y/ j (x,r+v J h) 

+ zJo exp[-Cj(r,x'H k (r+VjX',h - x')](v/v k ) 
x c. k ( r +Vj x') y r k [x, r+VjX' + v k (h-x')] dx' (38) 


To evaluate equation (38), we use the mean value 
theorem that guarantees linear terms of the final integral 
to be zero. First, we expand the attenuation factor as 

Q(r,x') = Io x ' Ojir+vj x") dx" 

« // [o/r+Vj h/2) + drO/r+Vj h/2) v } (x"- h/2)] dx" (39) 
and similarly for 

W+VjX',h - x') = f 0 h ~ x ' a k [r + vj x'+v k (h-x")] dx" 

~ So ~ x ' [Ojfr+Vj x'+v k h/2) 

+ d r a j (r+v J x'+v k h/2) v k (x"- h/2)] dx" (40) 

while applying the mean value theorem to the remaining 
factors of equation (38) and neglecting all but linear 
expansion terms in the integrand yields 

¥j(x+h,r) - exp[-Cj(r,h)] yz/xj+Vjh) 

+ Z k (v /vj a Jk (r+Vjh/2) y/ k [x, r+(Vj + v k )h/2] 
xf 0 h expfcjjfr+Vjhtyx'-crJr+fVj+VtJh&Xh-x')]} dx' 
- exp[-C/r,h)] y/j(x,r+Vjh) 

+ Z k (v/v k ) a jk ( r +Vjh/2) y/ k [x, r+(Vj + v k )h/2] 
x [exp{-Oj(r+Vj h/2)h}-exp{-a k [r+(v,+v k )h/2)]h}] 
/{(J k [r+(Vj+v k )h/2)] - (j/r+vj h/2)} + 0(h 2 ) (4 1 ) 

to be compared with the original HZETRN algorithm to 
0[(v j -v k )h] derived by Shinn et al. (1992) given as 

y/j(x+h,r) « exp[-^(r,h) ] y/j(x,r+Vjh)+I k (v/v k )cr Jk (r) 
xy r k (x,r+vjh) {exp[-Oj(r)h]-exp[-o k (r)h]} /]a k (r) - a/r)] (42) 

In earlier versions of BRYNTRN for proton/neutron 
transport, the flux scaling relation was taken correctly as 

yfj(x,r) = S(E) p/x, E) (43) 

but carried over to the final BRYNTRN for light- 
ions/neutron transport in the most recent version of 
BRYNTRN (Cucinotta 1993) 21 . In coupling to HZETRN 
with scaling given by 

y/(x,r) = V| S(E)p/x, E) (44) 

there is an inconsistency in flux scaling which must be 
accounted. The appropriate coupling is given in 
equations (38) through (42) with the added factor of vj/v k 
in the field coupling terms. The main effects on solution 
of the Boltzmann equation are expected for the light ions 
of H 2 , H 3 , and He 3 with only minor effects on the major 
light-ion/neutron components (n, H 1 , He 4 ). Note the main 
difference in the original and corrected code is 
determined by the ratios of vj/v k of the field coupling 
terms. 

To evaluate these differences in flux scaling, we have 
used the algorithm of equations (31) thru (33) for 
comparison with the original light-ion/neutron 
propagator, we use the 29 September 1989 solar particle 
event spectrum because of its relation to the 23 
February 1956 event (Wilson 2000) 22 represented by the 
proton spectrum (p/cm 2 -MeV) at the boundary 
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approximated above 30 MeV by (Nymmik 1997, 
Cleghorn and Badhwar 1999) 23,24 

<p p (0, E) = (2.034x1 0 7 /|3) x [p(E)/p(30)r 45 (45) 

where p(E) is the proton momentum (MV) given as 

p(E) - ?[E(E + 1876)] (46) 

and [3 is the proton speed relative to the speed of light. A 
low energy correction below 30 MeV mainly affecting 
transport results for depths less that 1 g/cm 2 in most 
materials is also added as 

0 (0, E) = 1.416 x1C?x exp[- p(E)/102. 1 18] 

x (E + 938)/p(E) (47) 

in good agreement with spectrometer data of the GOES 
satellite (Kim et al. 1999) 25 . The integral fluence values 
above 0.01 A MeV for neutrons, H 1 , and He 4 with Vj = 1 
are nearly unchanged (indistinguishable in Figs. 4-6) as 
they are the major components produced in reactions 
and H 1 is dominated by the fluence at the boundary over 
the first half of the mean free path. The H 2 and H 3 
integral fluences are decreased according to their Vj 
factors with values of 1/2 and 1/3 respectively (Figs. 7, 
8). The He 3 integral fluence is increased by the factor of 
Vj =4/3 as seen in Fig. 9. It is expected that dose will 
change little as the excess of doubly charged He 3 
contribution will largely cancel the singly charged H 2 and 
H 3 deficit contributions (approximately by a factor of (4/3 
- 7/6) times the total minor contributor’s dose). 



Fig. 4. Integral neutron fluence in aluminum shield using 
original and modified HZETRN code due to September 29, 
1989 solar particle event. 



Fig. 5. Integral proton fluence in aluminum shield using 
original and modified HZETRN code for September 29, 1 989 
solar particle event. 


The second correction to the propagator algorithm 
derived above, concerns the added accuracy of the HZE 
propagator to 0(1?) in equation (41) as opposed to the 
original HZETRN with error term 0[(Vj-v?)h]. We will 
show that the improved HZE propagator of 0(1?) allows 
control of the propagated error as well as reducing the 
local truncation error as we will demonstrate in the next 
section. 

Numerical Analysis of Marching Procedures 

There are two variables for which numerical 
approximation enter into the HZETRN propagator 
algorithms. The first is in the position variable xand the 
second is the residual range variable r. The coupling 
integrals of the Boltzmann equation involve integrals 
over energy that become principally integrals over 
residual range for the scaled flux equations although the 
energy shift operator of the Boltzmann equation couples 
residual range shift and position drift operators along the 
characteristic curves of the transport solution (Wilson 
1977) 11 . The principal concern is the necessary control 
of local truncation errors to insure that propagated error 
is controlled. In consideration of how errors are 
propagated, the error introduced locally by evaluation of 
y/(x, r ;■ + h) over the range (energy) grid with which it is 
defined as 

y/(x + h, rj = exp(- a h) \ff(x,r, + h) (48) 

whereas the local truncation error is given by 

W(x, r, + h) = y/Jx^ + h) + £,(h,r) (49) 

After the t? h step from the boundary, the numerical 
solution is 

V(kh, rj = exp(- a h) y J(k-1)h,r, + h] 

+ E X J' : exp[- o(k - X)h] £l (h,rJ (50) 

Suppose the local truncation error is bounded above 
such that e x (h,r J < e(h) for all X, then the propagated 
error is bounded by 

Sfjh) = E.J- 1 exp[- a(k - X)h] e x (h) 

< e(h) E^g' 1 exp[- o(k - X)h] ~ e(h) [1 -exp(-okh)]/ho (51 ) 

which is well behaved for all k and h if the local 
truncation error is bounded above by at least 0(1?) as 
demonstrated by Wilson et al. many years ago (1 991 ) 6 . 
The propagated error grows to a maximum of e(h)/ho 
requiring the 0(1?) limitation on the local error. The 
asymptotic bound for deep penetration is found to be 

t' vr p(h) < e(h) exp(-ah)/ [1- exp(-ah)] (52) 

emphasizing again the need to control the local 
truncation error as ho => 0. The original BRYNTRN and 
HZETRN propagator algorithms marginally met these 
requirements as was found by the detailed studies of 
Shinn et al. (1991, 1992) 20,26 . In the reductions leading 
to equations (31), (32) and (41), the error terms are 
0(1?) when the base algorithms are obtained but the 
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errors associated with the numerical approximation of 
the remaining functions of residual range (or energy) 
have been left so-far unspecified and were the subjects 
of prior studies (Shinn et al. 1991, 1992) 20,26 to be 
reviewed herein. 

The original BRYNTRN and HZETRN codes assumed 
approximate log-linear dependence of all discretized field 
functions of residual range that are on 0(A 2 ) for galactic 
cosmic ray like spectra where A is the order of the 
residual range spacing but only O(A) for most model 
solar particle events or trapped proton spectra (Wilson et 
al. 1 991 ) 6 . 

The original range-grid was derived using a uniform 
log(E)-grid of thirty points converted to range using 
range-energy relations of the transport media. Detailed 
studies by Shinn et al. (1 991 ) 26 used a 90-point log(E)- 
grid as standard for evaluation of errors in the original 
30-point grid and a 60-point grid. Maximum errors were 
first quantified to be a few percent in dose and dose 
equivalent at the largest depths of 150 g/cm 2 in air 
(Tables 2&3 of Shinn et al. 1 991 ) 26 . A systematic study 
of grid generation and numerical interpolation was then 
committed. It was found that a uniform log(r)-grid of 60- 
points gave an accurate interpolation (fraction of a 
percent of flux) with a fourth order Lagrange 
interpolation. It was desirable at that time to minimize 
the number of grid points as computational time is 
dominated by evaluation of the integral coupling terms 
and increases as N 2 . It was clear in the studies of 
uniform log-grid spacing that only the midrange errors 
were significant so that Shinn et al. replaced the fully 
uniform grid with a uniform grid over two sub-domains 
allowing even greater accuracy with only 30 grid points. 
An excess number of points are over the range of 1 
g/cm 2 with fewer points at the lower range values and is 
now a general feature of both BRYNTRN and HZETRN 
codes. The errors due to the residual range grid below 1 
g/cm 2 has no effect on the propagated error as the step 
size is on the order of 1 g/cm 2 so that this low energy 
part of the spectrum is deposited in the sub-range of the 
next step. This is facilitated by the scaled flux that 
approaches a constant at these lower energies [see 
equations (20) and (21)]. 

Aside from the issue of numerical interpolation and direct 
effects on the propagation routines is the evaluation of 
integrals of field quantities related to coupling terms. 
The original BRYNTRN and HZETRN codes used the 
assumed log-linear dependence and evaluated 
quantities analytically arriving at computationally efficient 
procedures (an important feature on contemporary 
machines at that time). Shinn et al. (1 991 ) 26 made 
detailed studies of numerical integration errors using the 
90-point solutions as a standard for which the original 
algorithms for integral flux resulted in errors of less than 
0.5 percent. It was found that substitution of a three- 
point Simpson’s rule reduced the integration errors by 



depth (g/cm 2 ) 

Fig. 6. Integral He 4 fluence in aluminum shield using original and 
modified HZETRN code for September 29, 1989 solar particle 

event. 



Fig. 7. Integral H 2 fluence in aluminum shield using original and 
modified HZETRN code for September 29, 1989 solar particle 

event. 



depth (g/cm 2 ) 

Fig. 8. Integral H 3 fluence in aluminum shield using original and 
modified HZETRN code for September 29, 1 989 solar particle 

event. 



Fig. 9. Integral He 3 fluence in aluminum shield using original and 
modified HZETRN code for September 29, 1989 solar particle 

event. 
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approximately an order of magnitude using midpoint 
values of the improved interpolation algorithm with the 
modified uniform log(r)-grid on two sub-domains. The 
reformulated propagation routines were found to have a 
fraction of percent error over the transport domain to 150 
g/cm 2 depths. In every case so far studied, the 
approximations in equation (41) are assumed correct 
and attention is given to evaluation of the right hand side 
without reference to the original integral on the left side 
of equation (32). The adequacy of the approximation in 
equation (33) is being addressed in ongoing activities. 

Shinn et al. (1 991 ) 26 examined the step size 
convergence within the BRYNTRN algorithm using the 
aforementioned modifications with the 30-point 
converged results. The step size was varied from 1 




(a) Aluminum shield on water. 20 g/cm^ Al ; 30 g/cm^ water. 

Fig. 10. Total dose and dose equivalent (ICRP 26) for the 
Webber benchmark SPE spectrum. 


g/cm 2 to 0.1 g/cm 2 for which dose for protons converged 
quickly but neutrons more slowly. The compromise step 
of 0.5 g/cm 2 is now standard in the BRYNTRN code and 
in the light ion propagator of HZETRN. The current 
version, so configured as discussed above with 30 
log(r)-grid points, results in 5 percent accuracy to 150 
g/cm 2 and is sufficient for most applications. Even so, 
standard practice now uses 80 such grid points assuring 
even improved accuracy for both GCR and SPE 
applications. Furthermore, the number of grid points is 
further adjusted to accommodate the simulation of 
geomagnetic cutoff effects while maintaining high 
numerical accuracy. 

In addition to the improved numerical procedures, it is 
noted that the light-ion/neutron propagator of the 1995 
HZETRN has substantial differences from its BRYNTRN 




(b) Iron shield on water. 20 g/crn^ Fe; 30 g/cm^ water. 

Fig. 10. Concluded. 


(version 3) predecessor. To quantify these differences, 
we evaluate dose and dose equivalent (as given by both 
the ICRP 26 and ICRP 60 quality factors) in 30 cm of 
water behind a 20 g/cm 2 shield of aluminum (and 
alternately iron) for the Webber (1963) 27 approximation 
of the 23 February 1956 spectrum (p/cm 2 -MeV) given as 
a P 0 ~ 100 MV spectrum with 10 9 protons /cm 2 above 30 
MeV in the following 

<!>J0, E) - 1(f x exp{[239. 1- p(E)]/100} 

x (2E + 1876)/[200 x p(E)] (53) 

and comparing with the Monte Carlo results (an early 
version of HETC with ICRP 26 quality factors) of Scott 
and Alsmiller (1968) 28 and more modern Monte Carlo 
codes using ICRP 60 quality factors. We first re- 
evaluate the 2005 HZETRN version with the ICRP 26 
quality factors by comparing to Scott and Alsmiller, the 
results are shown in Fig. 10 (a and b). A breakdown into 
specific contributions is presented by Wilson et al. (1991, 
1995) 6,29 using an earlier version of BRYNTRN (version 
1) and 1995 HZETRN, respectively. 

We have further approached code testing with a 
benchmark by neglecting the integral term of equation 
(32) and boundary condition given by equation (53) in 
both the analytic solution and 1995 HZETRN code. The 
analytical solution is given in equation (35) neglecting 
the integral term and unsealing the result according to 
equation (24). The initial testing of the 1995 HZETRN 
version chosen at random from various copies revealed 
that the light-ion/neutron cross section routines were 
corrupted and we replaced them by more accurate (and 
uncorrupted) routines developed by Tripathi et al. (1998, 
1999) 30,31 . Now the transported flux is generally within 1 
percent of the analytic solution as is the dose using 
Simpson’s rule but dose equivalent was found to be low 
by a few percent. Replacing Simpson’s rule by a ten- 
point Gauss-Legendre quadrature brings dose 
equivalent to within 0.15 percent of the analytic result 
and Gauss-Legendre quadrature will be a permanent 
feature of the 2005 HZETRN code with comparisons in 
Table 1. 

The percent differences of the analytical proton flux and 
the numerically generated proton flux at the iron- 
shield/water interface and at exit of the water slab are 
shown in Fig. 11. The results in Fig. 11 provide a direct 
test of the basic propagator methodology that is shown 
to be quite accurate. In addition to allowing evaluation 
of the accuracy of basic transport procedures and the 
nuclear attenuation factors, this benchmark provides a 
direct test of using equation (24) for unsealing the 
numerical solution developed on scaling relation (44) 
and demonstrating the requirements for the low energy 
equilibrium solution of equation (15) to be accurately 
maintained by the approximate numerical propagation 
method. The dose and dose equivalent at various 
depths in water for this analytic benchmark solution is 
given in Table 1 . This benchmark solution is now a code 


feature useful for validation after porting to other 
platforms and differing compilers. 

Table 1 . Dose and dose equivalent (ICRP 60) of penetrating 
protons from analytical solution compared to numerical 
solution (in parenthesis). Webber spectrum on 20 g/cm 2 of iron 
shielding 30 cm water. 


Depth, cm 

Dose, cGy 

Dose equivalent, cSv 

0 

8.405 (8.405) 

11.520 (11.505) 

5 

4.083 (4.074) 

5.009 (4.979) 

10 

2.321 (2.316) 

2.817 (2.800) 

15 

1.417 (1.414) 

1.707 (1.696) 

20 

0.909 (0.907) 

1.089 (1.082) 

25 

0.604 (0.603) 

0.720 (0.716) 

30 

0.412 (0.411) 

0.490 (0.487) 


* values in parenthesis are numerical solution 


A similar analytic benchmark has been developed for the 
1977 Solar minimum galactic cosmic ray spectrum and 
demonstrates that the propagator ignoring secondary 
particle production and fragmentation are a fraction of 
percent of the corresponding analytic solution with main 
errors near the boundaries of the energy grid as shown 
in Fig. 11 and most values are correct to a small fraction 
of 1 percent. The dose and dose equivalent of the 
analytic benchmark solution and numerical benchmark 
solution differ by less than 0.15 percent. 

Benchmarking can be important in both evaluation the 
code accuracy as well as provision of test cases for code 
verification after porting to other platforms and/or 
compilers. 



(a) analytic SPE benchmark 


%diff. in analytic vs. numerical solution (1977 GCR) 



(b) analytic GCR benchmark 

Fig. 1 1 . Numerical errors in proton spectra for analytic 
(a) SPE and (b) GCR benchmarks versus energy index. 
Indexed energies for SPE range from 0.01 to 900 MeV 
and for GCR from 0.01 to 50,000 MeV. 


The results of the BRYTRN (version 3), the 1995 
HZETRN (including ten years of drift), and the 2005 
HZETRN (improved numerical procedures as developed 
above) are shown in Fig. 12. There were many reasons 
for the differences including corruption of a nuclear 
reaction routine for light ions and a nuclear fragment- 
ation database in addition to development of improved 
numerical procedures. We have made appropriate 
modifications as discussed above resulting in the 2005 
HZETRN version with corrected nuclear routines and 
database. We have used a benchmark in the past 
based on an HETC result using the Webber spectrum for 
30-cm slab of water shielded by 20 g/cm 2 iron (or 
aluminum) as shown in Fig. 10 for the 2005 HZETRN in 
comparison with dose and dose equivalent (ICRP 26 
quality factor) according to HETC. It is clear that the 
result in Fig. 12 that considerable drift in codes have 
occurred over the years. 



depth in water, cm 


Fig. 12. Dose equivalent (ICRP 60) in water shielded 
from the Webber spectrum by 20 g/cm 2 of iron. 

For future purposes, the dose and dose equivalent in 
water are given in Table 2 for 20 g/cm 2 shields of 
aluminum and iron at the end of this paper. Values for 
the 1977 Solar minimum GCR spectrum for the 
aluminum or iron shielded water are shown in Table 3. 
In these tables, values for dose, expected TLD100 
response, and dose equivalent with ICRP 26 and ICRP 
60 quality factors are given. Any future code usage 
should reproduce these tables to ensure that the codes 
have not accumulated additional software errors, 
platform errors, or compiler related errors. 

Recently, new benchmarks have been provided for the 
two shield configurations described above (20 g/cm 2 of 
iron or aluminum shielding water) from the Monte Carlo 
Codes PHITS developed by the Japan Atomic Energy 
Agency (JAERI/JAEA) 32 and MULASSIS developed by 
the European Space Agency (ESA) 33 . The PHITS results 


are provided by Tatsuhiko Sato of JAEA and the 
MULASSIS results are provided by Jordi Bernabeu of 
Universitat Politecnica de Catalunya in collaboration with 
Hugh Evans of ESA/ESTEC. The Monte Carlo results 
for the Webber spectrum are shown in Table 4 to be 
compared with 2005 HZETRN results in Table 2. Dr. 
Sato also provided PHITS results for the 1977 Solar 
Minimum GCR spectrum given in Table 5. We note that 
differences between deterministic and Monte Carlo 
approaches tend to grow near the exit of the water 
column that may be caused by neutron (and lesser 
proton) leakage on the back surface that is not present 
in the current form of 2005 HZERN. There are other 
differences, especially for 1977 Solar Minimum GCR 
penetration problem, on the order of ten to twenty 
percent in dose and dose equivalent. 

Future HZETRN Development 

The next immediate step in code development is to 
validate this 2005 HZETRN using dosimetry and other 
instruments aboard shuttle and ISS for which improved 
dynamic/anisotropic environmental models have been 
recently derived (Wilson et al. 2005b) 34 In this activity, 
special emphasis will be given to validation and 
estimation of uncertainty similar to past activity (Wilson 
et al. 2005a) 5 Parallel to this validation activity is the 
development of the next generation engineering version 
of HZETRN with anticipated release in 2006. 

Recent research activity associated with HZETRN has 
had two parallel tracks. The first track is to advance 
Green’s function methods in anticipation of producing a 
code that is capable to be validated using high-energy 
ion beams (for example, Tweed et al. 2005). 35 The 
second track is to treat the off-axis scattering in the 
propagation of the light-ion/neutron propagator 
(Clowdsley et al. 2000, Heinbockel et al. 2003, Slaba et 
al 2006). 36,37,38 The next anticipated version of HZETRN 
will still use marching procedures for forward produced 
components of the interactions and evaluate the 
production source terms with broad angles with more 
appropriate angle dependent propagation techniques. 

CONCLUSION 

The reformulation of the HZETRN propagators with 
higher-order local truncation errors allows improved 
control of error propagation in the basic marching 
procedures. Conversion to dose and dose equivalent 
use improved numerical procedures based on a ten point 
Gauss-Legendre formulation. Analytical benchmarks 
are included in the 2005 HZETRN for code verification 
and in Table 1 of this report as a portable test. A 
benchmark with an early version of the HETC Monte 
Carlo code is contained herein in Fig. 10. A benchmark 
using the 2005 HZETRN code is given in Tables 2 and 3 
for future verification on various platforms and/or 
compilers and to guard against drift in future 


applications. Tables 4 and 5 contain new Monte Carlo 
benchmarks for evaluation of Tables 2 and 3 and for 
future reference. 
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Table 2. Dose (cGy) and dose equivalent (cSv) in a 30 cm water slab protected by aluminum or iron shield from the Webber solar 
particle event spectrum. 


Water 
Depth, cm 

Aluminum Shield 
Thickness of 20 g/cm 2 

Iron Shield 

Thickness of 20 g/cm 2 

X 

D(x), cGy* 

H(x), cSv** 

D(x), cGy* 

H(x), cSv** 

0 

7.09 (6.83) 

11.86 (11.56) 

9.18 (8.84) 

15.39 (15.12) 

5 

3.86 (3.75) 

6.06 (5.99) 

4.68 (4.54) 

7.32 (7.26) 

10 

2.36 (2.28) 

3.84 (3.75) 

2.77 (2.68) 

4.45 (4.37) 

15 

1 .53 (1 .48) 

2.53 (2.61) 

1.77 (1.71) 

2.95 (2.86) 

20 

1 .04 (1 .00) 

1.85 (1.79) 

1.18 (1.14) 

2.07 (1.99) 

25 

0.74 (0.71) 

1 .40 (1 .32) 

0.83 (0.78) 

1 .52 (1 .45) 

30 

0.54 (0.51) 

1.08 (1.02) 

0.60 (0.57) 

1.16 (1.09) 


‘values in parenthesis are expected for TLD100, “values in parenthesis are for ICRP26 quality factors 


Table 3. Annual dose (cGy) and dose equivalent (cSv) in a 30 cm water slab protected by aluminum or iron shield from the 1977 Solar 
Minimum GCR spectrum. 


Water 
Depth, cm 

Aluminum Shield 
Thickness of 20 g/cm 2 

Iron Shield 

Thickness of 20 g/cm 2 

X 

D(x), cGy* 

H(x), cSv** 

D(x), cGy* 

H(x), cSv** 

0 

20.9 (18.9) 

76.0 (66.8) 

22.0 (19.7) 

85.5 (75.7) 

5 

19.0 (17.5) 

58.2 (51.7) 

19.4 (17.8) 

64.9 (57.5) 

10 

18.3 (17.0) 

51.2 (45.8) 

18.6 (17.3) 

55.8 (49.8) 

15 

17.7 (16.6) 

46.5 (41.9) 

18.1 (16.8) 

49.9 (44.7) 

20 

17.3 (16.2) 

43.3 (41.8) 

17.6 (16.4) 

45.9 (41.3) 

25 

16.9 (15.9) 

41.1 (37.2) 

17.2 (16.1) 

43.1 (39.9) 

30 

16.5 (15.5) 

39.4 (35.7) 

16.8 (15.8) 

41.0 (37.1) 
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‘values in parenthesis are expected for TLD100, “values in parenthesis are for ICRP26 quality factors 


Table 4. Dose (cGy) and dose equivalent (cSv) in a 30-cm water slab protected by aluminum or iron shield from the Webber solar 
particle event spectrum evaluated using recent Monte Carlo codes: PHITS and MULASSIS/GEANT4 (within parenthesis). 


Water 
Depth, cm 

Aluminum Shield 
Thickness of 20 g/cm 2 

Iron Shield 

Thickness of 20 g/cm 2 

X 

D(x), cGy* 

H(x), cSv* 

D(x), cGy* 

H(x), cSv* 

0 

7.09 (6.82+1 .3%) 

10.9 (10.67+3.3%) 

9.21 (8.95+1 .2%) 

14.6 (14.12+2.8%) 

5 

3.90 (3.76+1 .8%) 

5.95 (5.62+4.8%) 

4.74 (4.54+1 .5%) 

7.16 (6.55+3.2%) 

10 

2.37 (2.27+2.2%) 

3.70 (3.48+7.2%) 

2.79 (2.72+2.0%) 

4.26 (4.14+6.5%) 

15 

1.53 (1.48+2.8%) 

2.44 (2.14+6.3%) 

1.76 (1.73+2.5%) 

2.74 (2.56+6.8%) 

20 

1.03 (1.02+3.4%) 

1.70 (1.62+8.3%) 

1.17 (1.15+3.2%) 

1.87 (1.80+8.9%) 

25 

0.717 (0.72+4.3%) 

1.21 (1.05+7.0%) 

0.806 (0.85+3.8%) 

1.32 (1.33+14.5%) 

30 

0.511 (0.51+5.3%) 

0.843 (0.87+18.3%) 

0.565 (0.60+4.8%) 

0.902 (0.94+9.9%) 


‘values in parenthesis are tentative GEANT4 results 


Table 5. Annual dose (cGy) and dose equivalent (cSv) in a 30 cm water slab protected by aluminum or iron shield from the 1977 Solar 
Minimum GCR spectrum evaluated using the recent Monte Carlo code: PHITS 


Water 

Aluminum Shield 

Iron Shield 

Depth, cm 

Thickness of 20 g/cm 2 

Thickness of 20 g/cm 2 

X 

D(x), cGy 

H(x), cSv 

D(x), cGy 

H(x), cSv 

0 

23.1 

69.9 

24.6 

83.9 

5 

22.0 

56.3 

22.5 

63.2 

10 

21.6 

49.2 

21.8 

53.3 

15 

21.2 

44.6 

21.3 

47.2 

20 

20.8 

41.1 

21.0 

43.1 

25 

20.3 

37.8 

20.4 

39.1 

30 

18.6 

32.6 

18.7 

33.5 
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